Now let's look at another class of computational problems, the so-called Stokes problem. These equations describe the flow of a liquid [1], [2].
First, let us try to understand the meaning of these equations. How can fluid flowing in a two-dimensional area be described? Suppose we are looking for a solution to the Stokes problem (liquid flow problem) in a square shaped area \( \Omega = [0,1]^2 \). At each point \( (x,y)\in \Omega \) of our area, we want to know the velocity of the liquid flow which is a vector and which assigns to each point two coordinates of the velocity field (velocities in the direction of the two axes of the coordinate system at this point) \( \Omega \ni (x,y) \rightarrow {\bf u(x,y) }=(u_1(x,y),u_2(x,y))\in R^2 \) (where bold \( {\bf u } \) stands for a vector of two components) and a scalar pressure field which assigns to each point one pressure value at that point \( \Omega \ni (x,y) \rightarrow p(x,y)\in R \).
So we have three unknowns at each point, two velocity components and one pressure component. Now, we need at least three equations to solve our problem. The Stokes equation in its general form is written as
\( - \Delta {\bf u } + \nabla p = f \).
Remembering that the velocity field has two components, and recalling what Laplasian and gradient mean, we get two equations in the so-called strong form
\( -\frac{\partial^2 u_1(x,y)}{\partial x^2} -\frac{\partial^2 u_1(x,y)}{\partial y^2} + \frac{\partial p(x,y)}{\partial x } = f_1(x,y) \)
\( -\frac{\partial^2 u_2(x,y)}{\partial x^2} -\frac{\partial^2 u_2(x,y)}{\partial y^2} + \frac{\partial p(x,y)}{\partial y } = f_2(x,y) \)
We get a third additional equation assuming the additional properties of the fluid for which we want to calculate the flow. Traditionally, the Stokes equation assumes that the fluid is incompressible and that there is no source or loss of fluid in the region. Mathematically, this condition is written in the form of the so-called divergence. We write specifically that the divergence is zero
\( \displaystyle{div} {\bf u}=0 \) where \( \displaystyle{div} {\bf u} = \nabla \cdot {\bf u} = \left( \frac{\partial }{\partial x}, \frac{\partial }{\partial y } \right) \cdot \left(u_1,u_2 \right) = \frac{\partial u_1(x,y)}{\partial x} +\frac{\partial u_2(x,y)}{\partial y } =0 \)
So we have the following system of equations to solve
\( -\frac{\partial^2 u_1}{\partial x^2} -\frac{\partial^2 u_1}{\partial y^2} + \frac{\partial p}{\partial x } = f_1 \)
\( -\frac{\partial^2 u_2}{\partial x^2} -\frac{\partial^2 u_2}{\partial y^2} + \frac{\partial p}{\partial y } = f_2 \)
\( \frac{\partial u_1}{\partial x} +\frac{\partial u_2}{\partial y } =0 \)
To solve the partial differential equation that is satisfied in a given area, we must additionally determine what is happening at the edge of that area. In this way, we will specify the problem that we want to solve. Any other boundary conditions give us a different problem and, consequently, a different solution. In our case, we introduce the so-called Dirichlet boundary conditions on the flow velocity, which tell us the values of the velocity field on the boundary of the region area.
Suppose our area is a bay on the bank of a river that flows uniformly from left to right, along the upper boundary of our area
\( \Omega \). We assume that the flow velocity on the upper edge of the area is
\( {\bf u}(x,y) = (1,0) \textrm{ dla } x \in (0,1), y=1 \).
For the remainder of the boundary, we assume that the flow velocity is zero
\( {\bf u}(x,y) = (0,0) \textrm{ dla } x \in \{0,1\}, y\in (0,1) \cup x\in(0,1),y=0 \).
We have certain boundary conditions defining the values of the velocity field on the boundary of the area. Such conditions are called Dirichlet boundary conditions. In the Stokes problem, we also calculate pressure, so we must make some assumption about pressure. We assume that in total the pressure over the entire area is equilibrated (the pressure sampled in different places by multiplication by the test function and integrated gives a total of zero)
\( \int_{\Omega} p(x,y) dxdy = 0 \).
We can now construct the so-called weak formulation of our problem, which will be suitable for simulation with the finite element method.
The weak formulation is obtained as follows. We integrate and multiply our equation by selected functions called test functions which we will use to average our equation in the area where these functions are defined. The main difference to the projection problem is that we now have a system of three equations, so we need three test functions, one for each equation, we denote them
\( v_1(x,y), v_2(x,y), q(x,y) \). We multiply the equations and integrate over the area
\( -\int_{\Omega} \frac{\partial^2 u_1(x,y)}{\partial x^2}v_1(x,y) dxdy -\int_{\Omega} \frac{\partial^2 u_1(x,y)}{\partial y^2}v_1(x,y) dxdy + \int_{\Omega} \frac{\partial p(x,y)}{\partial x }v_1(x,y) dxdy = \int_{\Omega} f_1(x,y) v_1(x,y) dxdy \)
\( -\int_{\Omega} \frac{\partial^2 u_2(x,y)}{\partial x^2}v_2(x,y) dxdy -\int_{\Omega} \frac{\partial^2 u_2(x,y)}{\partial y^2}v_2(x,y) dxdy + \int_{\Omega} \frac{\partial p(x,y)}{\partial y } v_2(x,y) dxdy= \int_{\Omega} f_2(x,y) v_2(x,y) dxdy \)
\( \int_{\Omega} \frac{\partial u_1(x,y)}{\partial x}q(x,y) +\int_{\Omega} \frac{\partial u_2(x,y)}{\partial y }q(x,y) dxdy =0 \)
In a similar way to the bitmap projection problem, each selection of a test function \( v(x,y) \) for averaging our problem in the area where the test function is defined gives us one equation. Various test function selections \( v(x,y) : v_1(x,y),v_2(x,y),q(x,y) \) will give us different equations
\( \int_{\Omega} \frac{\partial u_1(x,y)}{\partial x}\frac{\partial v_1(x,y)}{\partial x } dxdy +\int_{\Omega} \frac{\partial u_1(x,y)}{\partial y}\frac{\partial v_1(x,y)}{\partial y} dxdy +\int_{\partial \Omega} \frac{\partial u_1(x,y)}{\partial n} v_1(x,y) dS + \int_{\Omega} \frac{\partial p(x,y)}{\partial x }v_1(x,y) dxdy = \\ = \int_{\Omega} f_1(x,y) v_1(x,y) dxdy \)
\( \int_{\Omega} \frac{\partial u_2(x,y)}{\partial x}\frac{\partial v_2(x,y)}{\partial x } dxdy +\int_{\Omega} \frac{\partial u_2(x,y)}{\partial y}\frac{\partial v_2(x,y)}{\partial y} dxdy + \int_{\partial \Omega} \frac{\partial u_2(x,y)}{\partial n} v_2(x,y) dS + \int_{\Omega} \frac{\partial p(x,y)}{\partial y }v_1(x,y) dxdy = \\ = \int_{\Omega} f_2(x,y) v_2(x,y) dxdy \)
\( \int_{\Omega} \frac{\partial u_1(x,y)}{\partial x}q(x,y) +\int_{\Omega} \frac{\partial u_2(x,y)}{\partial y }q(x,y) dxdy =0 \)
where, after integration by parts, new terms appear - edge integrals calculated perpendicularly to the edge.
In our Stokes problem to calculate flow in a gulf adjacent to a river, we assume that forces
\( f_1(x,y)=f_2(x,y)=0 \). In the boundary condition, we assume that the speed on the upper boundary is
\( {\bf u}(x,y) = (1,0) \textrm{ dla } x \in (0,1), y=1 \). In other words, we introduce the Dirichlet boundary condition for the flow velocity here.
There is one technical problem with this formulation. Well, the flow field on the upper boundary is non-zero (directed to the right), while on the side edges (on the left and right boundary) the flow field is equal to zero. So, we have a jump from zero to one at the first coordinate of the velocity field. First of all, such a jump is not physical (it cannot be that the water at one point does not flow, and just next to it, it flows quite quickly, the transition from "not flowing" to "flowing" takes place in a smooth, regular manner). Secondly, if we would like a numerical solution to approximate using the finite element method, such a solution is presented as a linear combination of continuous polynomials. There are no such polynomials that jump from zero to one. Setting the problem this way will result in numerical singularities in the upper left and upper right edge of the area. To avoid this, in this case, the following trick is performed
\( g(x,y) = 1 \textrm{ dla } x \in (\epsilon, 1-\epsilon), y=1 \)
\( g(x,y) = \frac{x}{\epsilon} \textrm{ dla } x \in (0,\epsilon), y=1 \)
\( g(x,y) = \frac{(1-x)}{\epsilon} \textrm{ dla } x \in (1-\epsilon, 1), y=1 \)
\( {\bf u}(x,y)=(g(x,y),0) \textrm{ dla } x \in (0,1), y=1 \)
Then we actually have zero velocity in the corners of the area and there are no numerical singularities in the finite element calculations.
We then do the so-called "shift" of our problem. We divide the sought solution into two components
\( {\bf u}(x,y) = (h(x,y)+w(x,y),u_2(x,y)) \)
where \( h(x,y)=g(x,y)(1-y) \) is the so-called extension of our nonzero Dirichlet boundary to the whole area \( \Omega \).
So, we actually do our shift only in relation to the first component of the velocity field (because the second component is equal to zero). Then the first equation is changing
\( \int_{\Omega} \frac{\partial (h(x,y)+w(x,y))}{\partial x}\frac{\partial v_1(x,y)}{\partial x } dxdy +\int_{\Omega} \frac{\partial (h(x,y)+w(x,y))}{\partial y}\frac{\partial v_1(x,y)}{\partial y} dxdy +\int_{\partial \Omega} \frac{\partial (h(x,y)+w(x,y))}{\partial n} v_1(x,y) dS + \\ + \int_{\Omega} \frac{\partial p(x,y)}{\partial x }v_1(x,y) dxdy = 0 \)
A member with a given function \( h(x,y) \) is known and we can move it to the right-hand side \( \int_{\Omega} \frac{\partial w(x,y)}{\partial x}\frac{\partial v_1(x,y)}{\partial x } dxdy +\int_{\Omega} \frac{\partial w(x,y)}{\partial y}\frac{\partial v_1(x,y)}{\partial y} dxdy +\int_{\partial \Omega} \frac{\partial w(x,y)}{\partial n} v_1(x,y) dS + \int_{\Omega} \frac{\partial p(x,y)}{\partial x }v_1(x,y) dxdy = \\ = \int_{\Omega} \frac{\partial h(x,y)}{\partial x}\frac{\partial v_1(x,y)}{\partial x } dxdy +\int_{\Omega} \frac{\partial h(x,y)}{\partial y}\frac{\partial v_1(x,y)}{\partial y} dxdy +\int_{\partial \Omega} \frac{\partial h(x,y)}{\partial n } v_1(x,y) dS \)
We can calculate all the terms on the right-hand side. Another unknown \( w(x,y) \) entered on the left-hand side can be marked back as \( u_1(x,y) \); however, we must remember after solving our Stokes problem that we ought to add back the term with the boundary extension to it.
We introduce markings on our left-hand and right-hand sides
\( a({\bf u}, {\bf v}) =\int_{\Omega} \nabla {\bf u} : \nabla {\bf v}dxdy \)
where \( : \) means summing up all components
\( =\int_{\Omega} \frac{\partial u_1}{\partial x}\frac{\partial v_1}{\partial x}dxdy + \int_{\Omega} \frac{\partial u_1}{\partial y}\frac{\partial v_1}{\partial y}dxdy + \int_{\Omega} \frac{\partial u_2}{\partial x}\frac{\partial v_2}{\partial x}dxdy + \int_{\Omega} \frac{\partial u_2}{\partial y }\frac{\partial v_2}{\partial y }dxdy \)
\( b(p, {\bf v}) = \int_{\Omega} p div {\bf v} dxdy \)
\( = \int_{\Omega}p \frac{\partial v_1}{\partial x} dxdy+ \int_{\Omega} p \frac{\partial v_2}{\partial y } dxdy \)
\( f({\bf v} )= \int_{\Omega} \frac{\partial h(x,y)}{\partial x}\frac{\partial v_1(x,y)}{\partial x } dxdy +\int_{\Omega} \frac{\partial h(x,y)}{\partial y}\frac{\partial v_2(x,y)}{\partial y} dxdy +\int_{\partial \Omega} \frac{\partial h(x,y)}{\partial n } v_1(x,y) dS \)
Our problem can now be written as follows

\( b(q,{\bf u } ) = 0 \quad \forall q \)
Note that we have replaced a system of three equations with a system of two equations. Our two equations are "coded" in the first equation. Now, we can take the test functions \( {\bf v }=(v_1,0) \) to get the first solution, and \( {\bf v }=(0,v_2) \) to get the second equation.
Having thus described our problem, we can now assume that our sought solution
\( ({\bf u },p)= (u_1(x,y),u_2(x,y),p(x,y)) \) it has speed terms equal to zero on the boundary \( u_1(x,y)=u_2(x,y)=0\textrm{ dla }(x,y)\textrm{ na }\partial \Omega \). Likewise, we can assume that our test functions are also equal to zero on the boundary
\( v_1(x,y)=v_2(x,y)=0\textrm{ dla }(x,y)\textrm{ na }\partial \Omega \)
We now define the functions used for solution approximation and for testing.
We define what our group of elements will look like, and our B-spline basis functions spanning the elements, giving a knot vector along axis \( x \) and knot vector along \( y \). Here we refer to the third chapter, which describes the way of defining base functions with the help of knot vectors and B-spline functions.
Along the axis \( x \) we introduce a knot vector [0 0 0 1 2 3 4 4 4], likewise along the axis \( y \) we introduce a knot vector [0 0 0 1 2 3 4 4 4].
We obtained two bases of one-dimensional B-spline functions \( \{ B_{i,2}(x) \}_{i=1,...,6 } \) and \( \{B_{j,2}(y)\}_{j=1,...,6 } \). Then we will create tensor products from them \( B_{i,j;2}(x,y)=B_{i,2}(x)B_{j,2}(y),i,j=1,...,6 \).
Note that our area \( \Omega \) is stretched over a square \( [0,1]^2 \), similar to our 6 * 6 = 36 basis functions which results from the definition of the knot vector [0 0 0 1 2 2 3 4 4 4].
We now take our approximation of the velocity fields
\( u_1(x,y) \approx \sum_{i,j=1,...,6} u_1^{i,j } B_{i,2}(x)B_{j,2}(y) \)
\( u_2(x,y) \approx \sum_{i,j=1,...,6} u_2^{i,j } B_{i,2}(x)B_{j,2}(y) \)
and pressure
\( p(x,y) \approx \sum_{i,j=1,...,6} p^{i,j } B_{i,2}(x)B_{j,2}(y) \)
We also use combinations of linear B-spline functions for testing
\( v_1(x,y) \approx \sum_{k,l=1,...,6} v_1^{k,l } B_{k,2}(x)B_{l,2}(y) \)
\( v_2(x,y) \approx \sum_{k,l=1,...,6} v_2^{k,l } B_{k,2}(x)B_{l,2}(y) \)
and pressure
\( q(x,y) \approx \sum_{k,l=1,...,6} q^{k,l } B_{k,2}(x)B_{l,2}(y) \)
We select the basis functions \( ({\bf v },q) = (v_1,0,0)=(B_{k,2}(x)B_{l,2}(y),0,0) \textrm{ dla }k=l,...,6 \) and put into the equation ( 1 ). We get matrices
\( A1_{i,j=1,...,6;k,l=1,...,6} \int_{\Omega} \frac{B_{i,2}(x)}{\partial x}B_{j,2}(y)\frac{\partial B_{k,2}(x) }{\partial x}B_{l,2 }(y)dxdy + \int_{\Omega} B_{i,2}(x)\frac{\partial B_{j,2}(y)}{\partial y}B_{k,2}(x)\frac{\partial B_{l,2}(y)}{\partial y }dxdy \)
\( B1_{i,j=1,...,6;k,l=1,...,6} = \int_{\Omega} \frac{\partial B_{i,2}(x) }{\partial x } B_{j,2}(y) B_{k,2}(x)B_{l,2}(y) d,6) \)
Then we take base functions \( ({\bf v },q) = (0,v_1,0)=(0,B_{k,2}(x)B_{l,2}(y),0) \textrm{ dla }k=l,...,6 \) and put into the equation ( 1 ). We get matrices
\( A2_{i,j=1,...,6;k,l=1,...,6}=\int_{\Omega} \frac{\partial B_{i,2}(x)}{\partial x} B_{j,2}(y)\frac{\partial B_{k,2}(x)}{\partial x}B_{l,2}(y)dxdy + \int_{\Omega} B_{i,2}(x)\frac{\partial B_{j,2}(y)}{\partial y }B_{k,2}(y)\frac{\partial B_{l,2}(x)}{\partial y }dxdy \)
\( B2_{i,j=1,...,6;k,l=1,...,6} =\int_{\Omega} B_{i,2}(x)B_{j,2}(y) B_{k,2}(x)\frac{\partial B_{l,2}(y)}{\partial y } dxdy \)
Additionally, assuming base functions \( ({\bf v },q) = (0,0,q)=(0,0,B_{k,2}(x)B_{l,2}(y)) \textrm{ dla }k=l,...,6 \) and put into the equation ( 1 ). We get matrices
\( Q1_{i,j=1,...,6;k,l=1,...,6} =\int_{\Omega } \frac{\partial B_{i,2}(x)}{\partial x} B_{j,2}(y) B_{k,2}(x) B_{l,2}(y) dxdy \)
\( Q2_{i,j=1,...,6;k,l=1,...,6} =\int_{\Omega} B_{i,2}(x)\frac{\partial B_{j,2}(y)}{\partial y}B_{k,2 }(x)B_{l,2 }(y) dxdy \)
And the right-hand side
\( F1_{k,l=1,...,6}= \int_{\Omega} \frac{\partial h(x,y)}{\partial x}\frac{\partial B_{k,2}(x)}{\partial x } B_{l,2}(y)dxdy +\int_{\partial \Omega} \frac{\partial h(x,y)}{\partial n } B_{k,2}(x)B_{l,2}(y)dS \)
\( F2_{k,l=1,...,6}= \int_{\Omega} \frac{\partial h(x,y)}{\partial y}B_{k,2}(x)\frac{\partial B_{l,2}(y)}{\partial y} dxdy +\int_{\partial \Omega} \frac{\partial h(x,y)}{\partial n } B_{k,2}(x)B_{l,2}(y)dS \)
Our system of equations is as follows
\( \left( \begin{matrix} A1 & 0 & B1 \\ 0 & A2 & B2 \\ Q1 & Q2 & 0 \end{matrix} \right) \left( \begin{matrix} u_1 \\ u_2 \\ p \end{matrix} \right) =\left( \begin{matrix} F1 \\ F2 \\ 0 \end{matrix} \right) \)
To get an accurate solution of the Stokes equation requires additional stabilization of the finite element method. In the following modules, two methods of stabilization will be discussed, one based on the residual minimization method, the other based on the DG (Discontinuous Galerkin) method.